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Abstract 

The analytic continuation methods of complex scaling (CS), smooth exterior scaling (SES), and 
complex absorbing potential (CAP) are investigated for the supercritical quasimolecular ground 
state in the U^2+_Qf98+ gygtem at an internuclear separation of ii = 20 fm. Pade approximants 
to the complex-energy trajectories are used to perform an extrapolation of the resonance energies, 
which, thus, become independent of the respective stabilization parameter. Within the monopole 
approximation to the two-center potential is demonstrated that the extrapolated results from SES 
and CAP are consistent to a high degree of accuracy. Extrapolated CAP calculations are extended 
to include dipole and quadrupole terms of the potential for a large range of internuclear separations 
R. These terms cause a broadening of the widths at the o/oo level when the nuclei are almost in 
contact, and at the % level for R values where the ISu state enters the negative continuum. 



1 



I. INTRODUCTION 



The ability to produce super-strong electric fields in experiments, either by super-intense 
lasers or in collisions with highly ionized heayy atoms, opens the possibility to study nonlin- 
ear phenomena. It allows for detailed testing of the interface of relativistic quantum mechan- 
ics and quantum electrodynamics by probing electron dynamics at an energy scale where 
phenomena such as electron-positron pair creation may be detected against the background 
of more conventional atomic physics processes such as ionization and positron production 
by time-dependent fields [1]. 

An electron in a Coulomb potential represents a problem that has been studied in both 
relativistic and non-relativistic physics. As the potential depth increases with nuclear charge 
Z, the ground state's energy decreases as . Since the spectrum of the Dirac equation for 
a Coulomb potential is not bounded from below (due to the negative-energy continuum 
E < — meC^), a potential of sufficient strength can yield resonance states when the ground- 
state energy Eis < — mcC^. This causes the QED vacuum to become unstable and decay 
by pair-creation to a charged vacuum state [2]. The potential is called super-critical in this 
case. 

For a Coulomb potential, a single nucleus would need to be much more massive than any 
observed element {Z > 169 for finite-size nuclei), but when two large nuclei get close enough, 
the combined two-center potential can become supercritical [2, 3]. As the two nuclei get 
closer the energy of the quasi- molecular ground state (IScr) decreases and so the resonance 
energy is embedded deeper into the negative-energy continuum of the Dirac spectrum. These 
negative-energy electron states can be reinterpreted using CPT symmetry as positron states 
of positive energy. Therefore, by studying these supercritical resonances, the dynamics of 
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systems at the energy scale of particle creation can be explored. 

It is possible to build a matrix representation of the hermitian two-center Hamiltonian for 
a supercritical potential in order to measure the resonance parameters from the Lorentzian 
distributions in the energy spectrum or from the density of states [4] . More accurate results 
were obtained from analytic continuation methods which yield discrete resonance states with 
the following properties: (i) bounded, square- integrable eigenf unctions, i.e., elements of the 
Hilbert space of the original Hamiltonian; (ii) complex eigenenergies, E^, whose real and 
imaginary parts agree with the center and width of the Lorentzian shape in the spectrum 
of the original Hamiltonian. Note that the imaginary part is positive, i.e., = Ej-^g + ir/2 
for Dirac states with E^^s < — mgC^ [4]. 

Analytic continuation makes use of some parameter ( to turn the original Hamiltonian 
into a non-hermitian operator. This raises two problems: (i) one has to optimize ( to find 
the best approximation to the complex resonance eigenenergy; (ii) one worries about the 
effect of the unphysical ( parameter on the final result. The optimization in step (i) is 
performed by stabilizing some measure that depends on the eigenenergy as a function of ( 
(such as, e.g., the magnitude \Eii{()\). 

Recent work for the method of adding a complex absorbing potential (CAP) has demon- 
strated how to obtain more accurate results by extrapolating the complex £'/^(C)-trajectory 
to C = [5] . The present work extends this idea to the relativistic Dirac equation using the 
complementary methods of smooth exterior scaling (SES) and CAP. The stability and pa- 
rameter independence of the results is demonstrated. The extrapolation technique enables 
us to perform an extension of the calculations beyond the monopole approximation to the 
two-center potential. Results are given for the three coupled channels, k, = —1, 1, —2, for 
a large range of internuclear separation, i.e., the resonance width broadening due to dipole 
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and quadrupole interactions is calculated. 



II. THEORY 

A. Complex absorbing potential method 

The method of adding a complex absorbing potential (CAP) to the Hamiltonian has 
been used extensively in atomic and molecular physics [6, 7, 8, 9], and was put on firm 
mathematical grounds by Riss and Meyer [fO]. It was recently extended to the relativistic 
Dirac equation for Stark resonances [11], and for supercritical resonances [4]. It works by 
extending the physical Hamiltonian by an imaginary potential which makes the Hamiltonian 
non-hermitian. In the case of the Dirac Hamiltonian a CAP is added as a scalar giving, 

HcAv = H-tTim{r), (1) 

where 77 is a small non-negative parameter determining the strength of the CAP and (3 is 
the standard Dirac matrix [3]. The function W{r) determines the shape of the potential 
and is tailored for the specific problem. Currently, the most common use of a CAP is as 
a stabilization method: one solves the system on an equally spaced mesh of values of 77 
and takes the minimum of 
parameters [10]. 

In the present work we chose 



as the closest approximation to the true resonance 

»7=r;opt 



W{r) = Q{r -r^){r -r^Y (2) 

for the CAP, where 6 is the Heaviside function. The CAP parameter allows for the 
turn-on of the potential outside of the "bound" part of the wavefunction [12]. 
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B. Smooth exterior scaling method 

Smooth exterior scahng (SES) is an extension of the complex scahng (CS) method whose 
mathematical justification was developed, e.g., by W.P. Reinhardt [13], and Moiseyev [14]. 
CS introduces an analytical continuation of the Hamiltonian by scaling the reaction co- 
ordinate, r, by r — > re*^, and has been used extensively in atomic and molecular physics 
[15, 16, 17, 18]. CS was recently extended to the relativistic Dirac equation for supercritical 
resonances [4] and in 3D for a Coulomb potential with other short range potentials [19]. SES 
relies on the same justification as CS, but uses a general path in the complex plane that is 
continuous; when using non-continuous paths one refers to the method as exterior complex 
scaling (ECS) [14]. A simple path is obtained by rotating the reaction coordinate into the 
complex plane about some finite position, rg, instead of the origin. The transformation then 
has the form. 



It offers the advantage of turning on the scaling at a distance which can be chosen 
appropriately for a given potential shape. In analogy to the parameter of the CAP 
method it is natural to choose such that the "bound" part of the resonance state is not 
affected directly by the complex scaling. 

This additional freedom introduces some complications not found in CS: unlike CS, SES 
does not always have a minimum in the | | curve as a function of 6, which makes it difficult 
to determine a stabilized 9 value, ^opt, for Er{6). An approximation can always be made 
by finding the cusp of the trajectory in {i^rcsjE} space [20], although this does not allow 
for a very precise determination of ^opt compared with a minimum in | | . Alternatively, 
the parameters can be determined from either the or ^ curves, since in practice it 




(r — Tg) e*^ + for < r. 



r 



for r < Tg 



(3) 
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is usually found that at least one of them will have a minimum, with the best results for 
each parameter coming from its own derivative minimum [21]. Optimal results are obtained 
when the three derivatives have the same value of 6'opt yielding the same values for Ej-cs and 
r(lSa). 

Although CAP and SES appear as separate methods, they are not independent. It has 
been shown how a CAP can be transformed into CS [22] and that SES is related to CAP 
[23]. It has also been shown that a transformative CAP (TCAP) gives the same Hamiltonian 
as SES with the unsealed potential and a (small) correction term [24]. 

C. Pade approximant and extrapolation 

The goal of the analytic continuation methods is to make the resonance wavefunction, 
■j/^res, a bounded function by choosing an analytic continuation parameter ( > (^crit- Although 
z/^res is an eigenfunction of the physical, hermitian, Hamiltonian, it is not in the Hilbert space 
since it is exponentially divergent. For a sufficiently large critical value of C = Ccrit (called 
^crit for CS, and rjcrit for CAP), ipj-f.^ becomes a bounded function and is therefore in the 
Hilbert space of the physical Hamiltonian [14]. Taking the lim^^o always yields a real 

eigenvalue corresponding to H{( = 0)\I' = -E'_r(C = 0)\I^ since H{( = 0) is hermitian when 
acting on bounded functions. 

The authors of Ref. [5] proposed the following: instead of taking directly the C = limit, 
an extrapolation of a part of the trajectory, Eii{(), namely for ( > (cnt, is used to obtain 
the complex- valued -E'padG(C = 0)- The points used for the extrapolation are computed 
eigenvalues restricted to a region where ip^es remains in the Hilbert space (as represented by 
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the finite basis). For the extrapolation the Fade approximant is given by, 



^pad.(c) = . (4) 



where pt and qj are complex coefficients, and Np = 2{Ni + 1) is the number of points used 
in the approximant [25]. The extrapolated value of -EpadelC = 0) is given simply by po, and 
follows from a given set of (^-trajectory points, 

e, = {e, = + jAC), J = l...iVp} , (5) 

where d > Ccrit- For practical considerations about the added complications of the method 
(choice of e^, order of the extrapolation) we refer the reader to section IIIC. 



III. RESULTS 

To examine the smooth exterior scaling (SES) method and compare it to the simpler 
complex scaling (CS) method, and to the method of a complex absorbing potential (CAF), 
we look at the supercritical system of a single electron exposed to the field of a uranium {A = 
238) nucleus separated by 20 fm from a californium [A = 251) nucleus. The same system was 
explored in Ref. [4] using the CS and CAF methods and in Ref. [1] using phase-shift analysis 
and numerical integration. The nuclei are approximated as displaced homogeneously charged 
spheres with a separation of R between the charge centers. In the center-of-mass frame the 
monopole potential for each of the two nuclei with respective charge Zj, radius i?n^ and 
center-of-mass displacement -Rcm is given (in units of h = c = rric = 1, Z = Zi, = Rq^) 
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for r > r+ 

for r_ < r < r+ 

for r < r_ 
(6) 



where r± = Rcm ± -Rn [26]. The expression is obtained from the potential for a homoge- 
neously charged sphere displaced by Rqm along the 2-axis and then expanded in Legendre 
polynomials. 

To solve the Dirac equation for the potential given by Eq. (6), a matrix representation is 
constructed using the mapped Fourier grid method [27]. The radial coordinate is mapped 
to a new variable by, 

^ s0-4OOOarctan^ ^ 
(tt - 4>f 

were s is a scaling parameter which allows for the tailoring of the mesh points. The 
transformation maps r G [0, oo) to G [0, tt) which allows for a more efficient coverage of 
the relevant phase space. 



A. Pade extrapolation of SES trajectories 

In SES (as in the other methods) the energy Eii{9) forms a trajectory which depends 
on the calculation parameters {N,s,rs). Figure 1 shows the ^-trajectory for = 250, 
s = 400, = 0.01 and = 2 and rg = 3 respectively as a sequence of squares and circles 
respectively. The values are displayed as the deviation from a reference value given by 
Ercs = -1.757930695(8)mec2 and F = 4.084798(12)keV, which was obtained by a large-size 
CAP calculation using A^ = 3000, s = 400 and = 2. The uncertainties were determined 
from calculations with different basis parameters, s. 
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Figure 1: (Color Online) The r(lScr) resonance ^-trajectory from N = 250, rg = 2 (red squares) 
and Ts = 3 (blue circles) displayed as the difference from a reference calculation (cf. text) for 
the U-Cf system at i? = 20 fm in monopole approximation. The x's mark the resonance energy 



obtained by minimizing 



(19 



The A'p = 8 data points used for the extrapolations are bracketed 



by the solid symbols, and a spacing of A9 = 0.01 was used. 

In Fig. 1, the Efi{6) results for small 6 values start in the upper right hand and abruptly 
change direction after a few points. This happens in the vicinity of the reference value given 
by the point (0,0). As 6 increases beyond the cusp in the trajectory, E^{9) changes more 
slowly, and the ^-trajectory points become denser. Due to the finite number of collocation 
points A^, the actual value of 6'crit (the minimum 9 value such that ipj-^s is bounded and 
properly represented within the finite basis) is higher than theoretically predicted [14], and 
is close to the cusp-like behavior in the 6'-trajectory. 

Optimal, directly calculated, approximate values for the resonance energy are chosen by 
finding the minimum of I^^L^ [13, 14], and are displayed as black crosses for both 



opt 
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choices of Vg. The sohd hnes represent A^p = 8 Pade extrapolation curves (cf. Eq. (4)). The 
values of ei and e^p from Eq. (5), i.e., the bracketing points used to determine the Pade 
approximant, are indicated by filled symbols. The values of 6i = 0.28 for the Vg = 2 curve 
and 9i = 0.36 for the rg = 3 curve were chosen, because their -EpadolC = 0) points were 
found to be closest to each other. Using other starting values for the Pade curves would 
result in indistinguishable extrapolation curves on the scale of Fig. 1 (as long as the starting 
point would not be too close to the cusp in the 6'-trajectory). From the Pade curves one 
obtains the best approximation to the resonance parameters by selecting the end points 
(-E'padolC = 0)), located close to the the origin (0,0), i.e., close to the N = 3000 reference 
calculation result. 

Although very different trajectories from those shown in Fig. 1 were obtained for different 
Ts values, the Pade curves still extrapolated to nearly the same C = resonance energy, as 
long as Tg was chosen outside the "bound" part of the wavefunction and not too large, i.e. 
1.5 < Ts < 6. Larger values of rg pose difficulties for the finite- iV computation as it cannot 
span the large-r region properly. 

Simple complex scaling (CS) is equivalent to SES with rg = and, thus, falls outside of 
the range of acceptable rg values. It yields results with a much larger deviation from the 
reference result, even when Pade extrapolation is applied. 

The scaling parameter, s, introduced by the mapping of the radial coordinate in the 
mapped Fourier grid method in Eq. (7), plays a dominant role as far as numerical accuracy 
is concerned. In Fig. 2 the magnitude of the relative error in the width r(lS(T) is shown as a 
function of s for SES calculations with = 500, while using two reference CAP calculations. 
The latter were obtained from the stabilization method with a large basis (A^ = 3000) and 
yielded values of E'res = —1.757930695, P = 4.084798 for scaling parameter s = 400, and 
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Ercs = -1.757930702, T = 4.084808 for s = 750 respectively. Each SES data point was 
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Figure 2: (Color Online) Magnitude of relative error for Fade extrapolated widths AF(lScj) for 
N = 500 as a function of the scaling parameter, s, for the U-Cf system at i? = 20 fm in monopole 
approximation. Data shown are obtained using two reference calculations from the stabilized CAF 
method with N = 3000 and Tc = 2; red plus symbols: s = 400; blue squares: s = 750 (cf. text). 

obtained by finding the closest intersection of the Vg = 2,3,4 Fade curves from all starting 
6''s (for 6 > Ocrit), and taking the average of these three best values. The use of different Ts 
values within the acceptable range would cause changes too smah to see on this plot. The 
resonance position, i^rcs; is much more stable with respect to s and therefore not shown. 

The technique of combining information from different calculations in order to select 
the ideal subset of 6 points helps with the following problem. For each calculation with 
fixed value of we looked at extrapolated values for the resonance width as a function of 
the start value 6i, and observed the deviation from the average value. It was found that 
numerical noise present in the data was minimal for some range of 6i values (different for each 
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calculation). This noise was investigated as a function of extrapolation order A'p. Combining 
information from calculations with different values of Ts allowed us to eliminate two sources 
of random uncertainties associated with finite-precision input to a sensitive extrapolation 
calculation. 

The results demonstrate that the extrapolated {6 = 0) SES = 500 calculations exceed 
the precision of the stabilized (non-extrapolated) = 3000 CAP calculations. Table I gives 
the resonance parameter results (-Eres, T) for two different extrapolations, namely A'p = 4, 8 
for basis size A^ = 500, and mapping parameter 300 < s < 700, which corresponds to a 
subset of the data shown in Fig. 2. The standard deviation of the mean, 

i 

where x is the average x-value, are given separately for resonance position and width. The 
results are encouraging since extrapolation provides a marked improvement over the stabi- 
lized CS method discussed in Ref. [4]. 
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Table I: Averaged resonance position and width from the extrapolated SES method with basis size 
N = 500 using extrapolation orders Ap = 8 and A^p = 4. The values were averaged over the stable 
range 300 < s < 700. The standard deviation of the mean is given by a separately for position 
and width. 
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B. Pade extrapolation of CAP trajectories 

We have shown that a quadratic complex absorbing potential (cf. Eq. 2) gives better 
results than simple complex scaling (CS) when used as a stabilization method [4]. We find 
that the performance of SES represents an improvement over that of CS and it is competitive 
with CAP. 

Even for the same calculation parameters {N,s,rc), the ?7-trajectory for CAP is found 
to be different from the 6 trajectory in SES. Figure 3 shows a trajectory pair for identical 
calculation parameters, namely = 500, s = 600, and rg = Tc = 3 respectively. Even 
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-0.001 
-0.002 
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AE,es(10-6meC2) 

Figure 3: (Color Online) SES ^-trajectory with = 500, s = 600, = 3 (blue circles), and 
CAP ?7-trajectory for N = 500, s = 600, rc = 3 (red diamonds) displayed as the difference from a 
reference calculation for the IScr state in the U-Cf system at i? = 20 fm in monopole approximation. 
The X 's mark the approximate values obtained from the stabilization method. The reference value 
was obtained by Pade extrapolation of CAP calculations with N = 3000, Np = 4 for rc = 2,3,4 
yielding E.^s = -1.7579307062(5)mcc2 and F = 4.08481303(8)keV. 
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though the trajectories are rather different, the Pade curves extrapolate to rather close values 
(-E'padc('7 = 0)) as (i?Padc(^ = 0)) from the SES method, i.e. the extrapolated resonance 
parameters for both methods are very similar. When looking at plots of the N = 500 r] = 
extrapolation results for CAP as a function of s, as was done in Fig. 2 for SES, we observe 
almost identical results. In table II, the average and standard deviation of the mean, a 
(cf. Eq. 8), for the CAP results for N = 500, 300 < s < 700 are given for two different 
extrapolation sizes, A'p = 4, 8. 
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Table II: Same as in Table I, but for the extrapolated CAP method with basis size N = 500. 
C. Comparison of CAP and SES results 

Both methods (CAP and SES) show similar behavior with respect to different aspects of 
the calculation. Both are relatively insensitive to the number of Pade points, A^p, used in 
the extrapolation. It was found that the results were unchanged, to the relevant precision, 
for 4 < A^p < 12. The effects of changing the distance between the analytic continuation 
parameter (A^ for SES and Ar/ for CAP) were similarly found to be small compared with 
s. 

It was found that one could choose optimal start values for the Pade approximation 
by minimizing the deviation between Efi{C, = 0) from different or calculations. This is 
simple to implement, and works well for different values of A^p. In this way one obtains results 
that are independent of the starting point for the rotation (rg in SES) or the imaginary 
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potential {tc in CAP). The results are therefore very stable with respect to the analytic 
continuation parameters and depend only on the parameters from the mapped Fourier grid 
method. As shown in Fig. 2 for SES (a very similar graph was obtained for CAP), there is 
a range of s for which results are stable. For larger basis size the stable s-range increases 
making a judicious choice of s less important. 

We have averaged the results from both the SES and CAP methods over the stable s 
region, for a basis size of A^ = 500 in tables I & II. The standard deviation for the width (j(T) 
within either method is below 10~^, while the results differ at this level. It is, therefore, of 
interest to determine the reliability of the error estimate which is based upon basis parameter 
variations using larger- A^ extrapolated calculations. Comparison of the width results with 
such an estimate based upon A^ = 3000 CAP and SES calculations indicates that SES and 
CAP converge to the same value, closest to the CAP value given in table II. 

Concerning the most appropriate order for the extrapolations it is worth noting that the 
CAP Pade trajectories are rather straight in comparison with the ones for SES, and Ap = 4 
might be more appropriate in this case. Nevertheless, we find no systematic improvement 
when going to the lower-order approximation (which might be deemed more stable with 
respect to numerical noise in the trajectory points). 

IV. COUPLED-CHANNEL CALCULATIONS 

The ability to compute the resonance parameters to high precision with moderate basis 
size (e.g., A^ = 500) allows us to explore the effects of higher multipoles which are present 
in the two-center interaction. We are not aware of prior investigations of such coupled- 
channel resonance calculations of supercritical Dirac states. While the effect on the resonance 
position is expected to be small, the sensitivity of the width to computational details (cf. 
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the different results discussed in [4]) indicates that some broadening of the resonance may 



occur. 



To account for such two-center potential effects the wavefunction is expanded using spinor 
spherical harmonics, XK,fi, 



(9) 



which are labeled by the relativistic angular quantum number k, (analogous to / in non- 
relativistic quantum mechanics) and the magnetic quantum number /i [3]. The Dirac 
equation for the scaled radial functions, /(r) = rF{r) and g{r) = rG{r), then becomes 

{h = c = i), 

-^'^f^ = -(E-l)9.+ Y. (X.,.mr,/2)|X.>. , (10) 

K=±l 

A itoo 

+ = (i^ + 1)/.- E {x~.,,\V{r,R)\x--n,,)h , (11) 

K = ±l 

where V{r, R) is the potential for two uniformly charged spheres displaced along the z-axis, 
which is expanded into Legendre polynomials according to V{r, R) = Xl/^o ^K'^' R)Pi{cos9) 
[3]. The monopole term, Vq, is given explicitly in Eq. (6), and for the present work we include 
the coupling terms required for the k, = ±1, —2 channels (i.e., the Vi(r, R){x±K,fi \Pi\ X±r,h) 
dipole and V2(r, i?)(x±K,^j I-P2I X±R,tj) quadrupole terms). The k = 1, — 2 channels (-P1/2 and 
P3/2 respectively) have the strongest coupling to the k = —1 (S-states) and are therefore 
expected to have the largest impact on the supercritical ground state. 

Table 111 gives the ISa resonance parameters for different internuclear separations for 
both the one- and three-channel CAP calculations. The CAP calculations were performed 
using a basis size of = 500 per channel with 300 < s < 700. The CAP parameters of 
Tc = 2, 3, 4 were used to obtain the most stable 6 range for extrapolation, which was carried 
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Table III: Averaged IScr resonance position and width from the extrapolated CAP method with 
basis size N = 500, A^p = 4, for single-channel (monopole) and three-channel calculations as a 
function of separation, R, in the U-Cf system. The values were averaged over the stable range 
300 < s < 700 using Vc = 2, 3, 4. The value in parenthesis represents the standard deviation from 
the average and the last column gives the relative difference of the width from three- to one-channel. 

out with order A^p = 4. Values in parentheses indicate the standard deviation of the mean 
(cf. Eq. 8). We note that the accuracy of the calculations - as indicated by the deviation - 
is much higher than what is required to measure the effect of the the P-state channels on 
the ISo" resonance. 
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The final column contains the relative difference of the width between the three-channel 
results and the monopole approximation. The correction due to the dipole and quadrupole 
potentials is seen to increase with internuclear separation R. This trend can be expected, 
as the overall interaction becomes less spherically symmetric in this limit. The effect is 
small, however, as the strongest contribution towards the binding energy (and thus the 
supercriticality) is provided by the monopole part. The dipole interaction contributions 
to the resonance parameters are a result of the relatively small charge asymmetry in the 
U92+_Qf98+ potential. Quadrupole couplings to D states can be expected to generate more 
significant changes in E^^s and T. 

V. DISCUSSION 

We have extended the method of smooth exterior scaling (SES) to the relativistic Dirac 
equation. For both the addition of a complex absorbing potential (CAP) and SES, we have 
also shown that by extrapolation of the complex energy eigenvalue trajectories which are 
functions of the analytic continuation parameter, (, to ( = using appropriate ^-trajectory 
values (cf. section II C), one obtains an accurate and robust estimation of the resonance 
parameters. These results use much smaller basis than would be required for similar pre- 
cision using stabilization. While adding the extrapolation adds new parameters, we found 
minimizing the distance between extrapolated results for different calculation parameters 
was a simple and highly effective method of reducing the parameter space. The effects of 
the number of calculation points, A'^p, was found to be very small provided it was in the 
range 4 < A'p < 12. 

The effects of higher order terms and the P-states were explored and found to have a in- 
creasing influence as the two-center system becomes less spherically symmetric which occurs 
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at larger internuclear separations R. The resonances could acquire additional broadenings 
and shifts from dynamical effects not included in the quasi-static approximation. 
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